The T = RFIM on a Bethe lattice: correlation functions along the hysteresis loop 
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We consider the Gaussian random field Ising model (RFIM) on the Bethe lattice at zero temper- 
ature in the presence of a uniform external field and derive the exact expressions of the two-point 
spin-spin and spin-random field correlation functions along the saturation hysteresis loop. To com- 
plete the analytical description and suggest possible approximations for the RFIM on Euclidian 
lattices we also compute the corresponding direct correlation functions (or proper vertices) and 
show that they decay rapidly with the distance in the weak-coupling/large disorder regime; their 
range, however, is not limited to the nearest-neighbor distance. 
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I. INTRODUCTION 

The random field Ising model (RFIM) at zero tempera- 
ture is a simple prototype of a class of disordered systems 
(such as random magnets, martensitic materials, fluids 
in porous solids,...) that exhibit hysteretic and jerky be- 
havior when slowly driven by an external field lj. The 
most interesting feature of the model which has been the 
subject of extensive analytical and numerical studies is 
the existence of a disorder-induced nonequilibrium phase 
transition between two different regimes of avalanches [2]. 
This transition manifests through a change in the shape 
of the magnetization hysteresis loop that evolves from 
continuous to discontinuous as the disorder strength is re- 
duced. The discontinuity is associated to a macroscopic 
avalanche involving a finite fraction of the spins in the 
thermodynamic limit. As shown recently, this type of 
mechanism plausibly explains the hysteretic behavior of 
4 He adsorbed in high porosity silica aerogels [3JE]. Inter- 
estingly, this nontrivial behavior is already present on the 
Bethe lattice (i.e. the infinite Cay ley tree) where a fully 
analytical characterization of the major and minor hys- 
teresis loops, the avalanche size distribution, and other 
quantities can be obtained thanks to the tree topology 5 - 
9 . In this case, the out-of-equilibrium phase transition 
occurs when the coordination number z > 4 and is de- 
scribed by a traditional saddle-node transition in the self- 
consistent field equation [5 (which makes the critical be- 
havior the same as that for the infinite-range mean-field 
model). In this work we extend the analytical descrip- 
tion to the spin-spin and spin-random field correlation 
(or Green's) functions along the hysteresis loop, using 
the fact that correlations on a tree-like graph have a one- 
dimensional character. For a Gaussian distribution of the 
random fields, the spin-random field correlation function 
is also related to the slope of the magnetization curve 
through a 'susceptibility' sum-rule. The motivation for 



this calculation is twofold. First, on the theoretical side, 
Green's functions (or, better, their matrix inverse, the so- 
called direct correlation functions in liquid state theory or 
proper vertices in field-theoretic language) may be used 
as the building blocks of approximate theories, as illus- 
trated by the recent computation of the hysteresis loop in 
the three-dimensional soft-spin random field model [TU]. 
Exact results, even for simple models, may give some in- 
sight of the actual structure of these functions. Secondly, 
on the experimental side, scattering methods are now fre- 
quently combined with other standard probes (response 
to an applied field or thermodynamic measurements) for 
extracting information on the structure and the dynamics 
of systems with quenched randomness (see e.g. Ref. |llj 
in the case of fluids adsorbed in porous solids) . Knowing 
the structure of the correlation functions can thus make 
easier the interpretation of the scattered intensity [12]. 

The outline of the paper is as follows. In section II, we 
define the model and give the expressions of the corre- 
lation functions, first for the one-dimensional chain (cor- 
recting the result obtained in Ref.[T5]). and then gener- 
alizing to the Bethe lattice (the detailed calculations are 
presented in Appendices A and B). Analytical predic- 
tions are compared to simulations performed on regular 
random graphs. In section III, we compute the corre- 
sponding direct correlation functions. We then conclude. 



II. MODEL AND CORRELATION FUNCTIONS 



The RFIM is defined by the following Hamiltonian 



H = -JJ2 S i S i + J2( H + W 



(1) 
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where the N spins Si — ±1 are placed on the vertices 
of a Bethe lattice with coordination number z. The first 
sum is restricted to nearest-neighbors (n.n.) pairs and 
J > 0. H is a uniform external field and the fields {hi} 
are random variables drawn independently from a Gaus- 
sian distribution p(h) — exp(— /i 2 /2A)/v2~7rA with the 
variance A measuring the strength of disorder. 



The relaxation dynamics is the T = limit of the 
Glauber dynamics and consists in aligning the spins with 
their local effective field at each time step[3], 



Si = sgn(/j) 



where 



dH 
dS, 






H 



(2) 



(3) 



and the sum runs over the z nearest neighbors of site 
i. The dynamics thus proceeds via a series of avalanches 
which stop when a metastable state is reached, i.e., when 
all spins satisfy Eq. pi. The saturation hysteresis loop 
is obtained by adiabatically ramping H from — oo to +oo 
and back. Thanks to the tree topology of the Bethe lat- 
tice and the abelian property of the dynamics (e.g. the 
fact that the metastable state after an avalanche does not 
depend on the order in which the spins flip), the shape 
of the hysteresis loop can be exactly derived. According 
to Ref.[5], the magnetization m(H) along the lower half 
(ascending) branch is given by 



1 



fc=0 



m(H) + 1] = E (l)p*(H) k [l - P*(H)]*- k p k (H) 



(4) 

where Pk(H) (k = 0..z) is the probability for a down spin 
to flip up at the field H when k of its z nearest neighbors 
are up, 



Pk(H) = 



+oo 



(z-2k)J-H 



p(h)dh = -erfc 



(z-2k)J-H 



(5) 

(here erfc(x) = (2/y/n) j°° duexp(— u 2 ) is the comple- 
mentary error function), and P*(H) is solution of the 
self-consistent equation 



P*(H) = E 



k=0 



Z-\ 

k 



P* (H) k [l - P* (H) 



-l-k 



Pk(H) 



(6) 

This key quantity represents the conditional probability 
that a nearest neighbor of spin i flips up before spin i. 
For z > 4, the polynomial equation p| has several so- 
lutions at low enough disorder (for A < A c (z)) and the 
magnetization displays a jump discontinuity at a coercive 
field H C (A). 

In the following we are interested in calculating the 
correlations along the loop between the spin at site i and 
the spin or the random field at site j, 



GSS Q Q Q Q 
ij — OiDj Ji Jj 

Or,-,,- = Jjtlj 



(7) 



where the overbar denotes the average over the random 
field distribution p(h) and the dependence on the applied 
field H is implicit (hence Si = m(H) as given by Eq.|4J) 
along the ascending branch). Due to the average over 



disorder, the two functions only depend on the distance 
between the two spins, i.e. on n, the number of bonds 
between i and j. We thus denote them by G 3S (n) and 
G sh (n), respectively. 

Let us recall that at finite temperature and equilib- 
rium, because of the additional average over thermal 
fluctuations, there are two distinct spin-spin correlation 
functions, < SiSj > — < Si >< Sj > and < SiSj > — 
< Si > < Sj > where < ... > denotes the thermal 
average[H]. The former (the so-called connected or trun- 
cated function) may be non-zero at T = if the ground 
state of the system is highly degenerate. This does not 
occur when the random-field distribution is continuous 
and then only the disconnected function G ss (n) remains 
non-zero. At T = 0, one may also consider an aver- 
age over all the metastable states at a given field H 
and then distinguish again connected and disconnected 
contributions [ID] . However, the connected contribution 
vanishes along the hysteresis loop since there is only one 
metastable state and, again, only G ss (n) remains. On 
a regular Euclidian lattice, its Fourier transform is the 
structure factor S(q) which is the quantity measured in 
scattering experiments. G ss (n) should not be confused 
with the avalanche correlation function that measures the 
probability that the initial spin of an avalanche will trig- 
ger, in the same avalanche, another spin a distance n 
away [I]. In particular, in finite dimension, the algebraic 
decays of these two functions at criticality are not de- 
scribed by the same exponent. 

As was noticed only recently [TU], for a Gaussian distri- 
bution of the random fields, there exists a 'susceptibility 
sum- rule' that relates the correlation function G sh (n) to 
the slope of the magnetization curve at T — 0. It is ob- 
tained by using the following property of the Gaussian 
distribution: 



dhp{h)hA(h) 



A / dh^A(h) 
8A(h) 



A / dhp(h) 



dh 



Hence 



bitlj 



dhj ' 



and by summing over i and j one gets 

-Vg s,1 = a ( 



'•3 



dm 
'dH 



(8) 



(9) 



(10) 



On the Bethe lattice, this becomes 



dm 
'dH 



(11) 



G sh (0) + J2cnG sh (n) = A- 

n=l 

where c„ = z(z — l) n_1 is the number of sites distant 
from an arbitrary site i by n > 1 bonds (i.e. the number 
of sites that belong to nth shell) . 



To compute the correlation functions we first consider 
the case of a ID chain (i.e. z = 2) and then extends 
the results to the Bethe lattice with generic coordination 
number z. We find that 



G ss (n) = X n - 1 [a + b(n-l) 
G sh {n) = \ n - l G sh {l) 



(12a) 
(12b) 



for n > 1 (with G ss (0) = 1 — m 2 due to th e ha rd-spin 
condition Sf = 1, and G sh (0) given by Eq. Q). The 
explicit expressions of A, G ss (l) = a, G ss ( 2) = X(a + 
b), an d G s/l (l) are given by Eqs. ([23|, JBlf , | |B10| and 
(B14), respectively. 



One dimension 




FIG. 1: (Color on line) Correlation functions G 33 (n) and 
G ah (n) along the ascending branch of the hysteresis loop for 
z — 2 and A = 4. The simulation results (symbols) are 
compared to the predictions of Eqs. (12) (lines). Simulations 
were performed on random graphs with N = 10 6 and the 
results were averaged over 1000 disorder realizations. 



The hysteresis loop in the ID chain was calculated in 
Ref.|15j. In this case, there are only three probabilities 
Po,Pi,P2 defined by Eq. ^, and from Eqs. Q and (pi) 
the magnetization along the ascending branch is simply 
given by 



2 [(1 - P*) 2 p + 2P*(1 - P*) Pl + P* 2 p 2 ] 



1 
(13) 



with P* = po/(l — p\ + po) (hereafter, to simplify the 
notation, the dependence of all quantities on the field H 
is dropped) . The results on the descending branch can be 
obtained by symmetry. The analytical calculation of the 
spin-spin correlation function G ss (n) was first considered 
i n R ef. [13] but the final expressions of a and b in Eq. 
(12a I are wrong (see Appendix A); the dependence on 
the distance n, however, is correctly given. Moreover, 



G sh (n) was not considered. On the other hand, the exact 
expressions of G ss (l) = SiS i+1 - m 2 and G sh (0) = Wh 
were derived in Ref. |7] in order to compute the energy per 
spin along the hysteresis loop. The complete calculation 
that leads to Eqs. (12) is performed in Appendix A. In 
particular, we obtain 



A = p x - po 



(14) 



as correctly found in Ref. |13j . As shown in Fig. 1, Eqs. 
(12) are in excellent agreement with the results of nu- 
merical simulations performed on random graphs with 
z = 2 (we also checked numerically that the expression 
of G ss (n) given in Ref.jT3~j is not valid). 

The two functions G ss (n) and G sh (n) are thus charac- 
terized by the same correlation length S,\ = (— In A) -1 = 
[— \n(px — po)] -1 - However, G ss (n) is not a purely expo- 
nential function because of the prefactor b(n — 1). Re- 
markably, a similar behavior has been observed for the 
equilibrium RFIM in the very few cases where the corre- 
lation function G s e s q (n) — < S l S i ± n > — < Si > < Si± n > 
has been calculated exactly. This is indeed the leading 
long-distance behavior observed at T > with the special 
random-field distribution (somewhat related to percola- 
tion) considered in Ref. [16] (in this model, however, the 
T = behavior is complicated and the correlation func- 
tion behaves at long distance as an exponential divided 
by n 2 [l?j). For the Gaussian distribution that we here 
consider, no analytical expression is available for generic 
values of A and T, but an exact result has been ob- 
tained in the universal regime where the random field 
and the temperature are both much smaller than the 
exchange coupling [TB]. For H = 0, the leading long- 
distance behavior turns out to be also proportional to 
nexp(-n/£ eg ), where £, eq = 8J 2 /(tt 2 A) = (2/ir 2 )L IM 
and Ljm is the Imry-Ma length that sets the typical size 
of the domains in the ID chain at T = 0[19]. This coinci- 
dence is noteworthy but it must emphasized that the full 
expression of G s eq (n) in this regime is much more com- 
plicated than the one described by Eq. (12a) (moreover, 



for H =/= 0, G s eq (n) decays as a sum of exponentials). The 
correlation length £>, along the hysteresis loop also be- 
haves quite differently from £ e? in the limit A« J: it 
goes to the finite value l/ln(2) in zero applied field (as 
H = does not play any special role along the hystere- 
sis loop) and grows like £\ ~ (^/ttJ/V^A) exp( J 2 /2A) 
for H = J, which is the value of the field for which the 
susceptibility dm/dH is maximum. 



An interesting consequence of Eq. (12a) is that the 



structure factor S(q) in the small-g regime is a superpo- 
sition of a Lorentzian and a Lorentzian-squared terms. 
By definition 



S(q) = G BS (0) 



OO 



}G ss (n) 



(15) 



and using 



+00 



£V9'Al'l = T 



xcosq 



with 



2A 



1 + A 2 



we obtain after simple algebra 

B 



S(q) =A + 



C 



1 — xcosq [1 — xcosq] 2 



(16) 



(17) 



(18) 



with 



A = l- 



D 



(a - b)Vl - x 2 



A 



C = b 



1-x 2 
A 



(19) 



The Lorentzian plus Lorentzian-squared structure that 



emerges from Eq. ( 18 1 in the small-q regime is also found 
in the mean- field theory of the equilibrium RFIM[5D] and 
is usually used to fit experimental data on random mag- 
netic systems. 

In contrast, the spin-random field correlation function 
G sh (n) is a pure exponential for n > 1 so that its Fourier 
transform simply reads 



G sh {q) = [G s "(0) - 



G s ' l (l). G sh (l) Vl-x 2 



A 1 — x cos q 



(20) 



This yields 



G s 



l (q = 0) = G sh (0) + - -G sA (l) , (21) 

1 — A 



and using t he e xpress ion o f the magnetization, Eq. Q, 



and Eqs. (B3| and (B14| for G sh (0) and G sh (l), one 



can check that the susceptibility sum-rule, G sh (q = 
0) = A(dm/dH), is indeed satisfied. Note that the q- 
independent term inside brackets in Eq. (20) is non-zero, 



which is a somewhat unusual feature (the constant term 
A in Eq. (18) is also non-zero because Eq. (12a) is only 



valid for n > 1). As will be discussed in more detail in 
section III in the case of the Bethe lattice, this has a sig- 
nificant consequence for the matrix inverse of G sh (the 
so-called direct correlation function). 



B. Bethe lattice 

In principle, the probabilistic reasoning used in Ap- 
pendix A for the one-dimensional chain can be extended 




FIG. 2: (Color on line ) Correlation functions G 33 (n) and 
G (n) on a Bethe lattice with coordination number z — 3 for 
A = 9 (the curves result from an average over 5000 random 
graphs of size N = 10 5 ). The simulation results (symbols) 
are compared to the predictions of Eqs. (12) (lines). 
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FIG. 3: (Color on line) Same as Fig. [5] for z = 4 and A = 4. 



to the case of the Bethe lattice with generic coordina- 
tion number z. This is how the analytical expressions 
of G ss (l) and G sh (0) were derived in Ref.[S] in order to 
compute the energy per spin along the hysteresis loop 
(thereby generalizing the ID results of Ref.[7]). These 
expressions are recalled in Appendix B where we also 
calculate G ss (2) and G sh (l). Bowever, using the same 
method to derive the general expressions of G ss (n) and 
G sh (n) is unnecessarily complicated. Instead, one can 
simply exploit the fact that there is a unique path con- 
necting a given pair of spins on a Bethe lattice so that 
the dependence of the correlation functions on the dis- 
tance n must be the same as in one-dimension (just like 
in nonrandom systems). This implies that Eqs. (12) are 
also valid for the Bethe lattice. Strictly speaking, we do 



not provide a demonstration of this assertion 1 but it is 
fully supported by numerical simulations for small values 
of n, as illustrated in Figs. 2 and 3. 

The analytical expression of A can then be obtained 
via the susceptibility sum-rule. Inserting Eq. (12b) in 
Eq. (fill yields 



A S = ="■<»> + T3V,,. 



— G">(1) (22) 



so that 2 



A = 



z-V 



\l-z 



G 8h (l) 



Adm/dH -G sh (0)> 



(23) 



Using Eqs. Q, ( |B3| and ( |B14| , one can check that Eq. 
( 14 ) is recovered for z = 2. 
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FIG. 4: (Color on line) The quantities \(H) (top) and n(H) 
(bottom) characterizing the exponential decay of the correla- 
tion and direct correlation functions, respectively, for z — 4 
and (a) A = 2, (b) A = A c w 3.173, and (c) A = 6. In (a), 
below the critical disorder A c , A and fi jump discontinuously 
at the coercive field H C (A). In (b), A -> (z - 1) _1 = 1/3 
at the critical field H C (A C ) = J. Note that /i is significantly 
smaller than A in (c), above the critical disorder. 



For z > 4 and A < A c (z), the magnetization jumps 
discontinuously at the coercive field H C (A), which cor- 
responds to a spinodal singularity where dm/dH — > 
+00 EJ. As can be deduced from Eq. ( [22] ) , this is 
due to the fact that A — > (z — l) -1 , as illustrated in 



Eq. I ]12a| for G ss (n) is actually confirmed by the very recent 
analytical calculations of Ref. [21|. In that work, however, there 
is some confusion between G BS (n) and the so-called 'avalanche 
correlation function' (using the terminology of Ref.fl]). This 
latter function can be shown to behave as a simple exponential 
for n > 1, without the n — 1 prefa ctor. 

Using Eqs. J4J, ( |B3| ) , and l |B14] , it can be checked that Eq. 
( |_'i| is equivalent to the compact expression obtained in Ref. l21| : 
A = (z - l)- 1 8F(P*)/dP*, where F(P*) is the the r.h.s. of Eq. 

15). 



Fig. 4b for z = 4. Therefore the correlation length 
£a = (— InA) -1 keeps a finite value at all H and A, in- 
cluding at the critical point. This is indeed the expected 
(and standard) behavior on the Bethe lattice where the 
divergence of the susceptibility is generated by the expo- 
nential growth of the number of sites at the distance n 
(due to the hyperbolic-like geometry of the lattice) and is 
not associated to a divergence of the correlation length. 

III. DIRECT CORRELATION FUNCTIONS 



We now investigate the structure of the direct cor- 
relation functions (or one-particle irreducible functions, 
or else proper vertices in field-theoretic language) which 
(roughly speaking) are the matrix inverse s o f the cor- 
relation functions Gfj 1 and Gfj 1 (see Eqs. Ety and (J46J) 
below). As briefly mentioned in the introduction, the mo- 
tivation for this calculation is that proper vertices may be 
simpler or at least shorter-ranged than the Green's func- 
tions, and therefore can be used as the building blocks of 
approximate theories. For instance, in liquid-state the- 
ory, the direct correlation function c{r) , which is the ma- 
trix inverse of the pair correlation function h(r) at equi- 
librium and is defined via the so-called Ornstein-Zernike 
(OZ) equation, is in general shorter-ranged than h(r) 
(essentially having the range of the pair potential), ir- 
respective of the thermodynamic state [52]. This feature 
is the starting point of the successful integral equation 
approach to the structure and thermodynamics of sim- 
ple liquids. For Ising spins on a lattice with nearest- 
neighbor interactions, it is also a reasonable approxima- 
tion to assume that the matrix inverse of the spin-spin 
correlation function is zero for n > l|23j. This can been 
used to build a very accurate description of the three- 
dimensional Ising model[23] and other spin models 25 , 
including in the presence of quenched disorder [26- More 
recently, a similar approximation has been proposed to 
obtain an analytical description of the hysteresis loop 
in the three-dimensional soft-spin random field model at 
T = [10] . It is therefore interesting to check whether 
the direct correlation functions on the Bethe lattice are 
indeed shorter-ranged than G sh and G sh . 

We first consider the function C sh — {C^} defined by 



the OZ equation ^2 k C ik 
C sh 



sh/^ish 



i.e. 



\G 



sh] 



(24) 



where G sh and C sh are N x N matrices. As pointed 
out in the preceding section, G sh (n) is a pure exponen- 
tial for n > 1, but G sh {l) ^ \G sh (0). At first sight this 
is an innocuous feature but it has an important conse- 
quence for C sh (n). Indeed, as is well known (and is also 
shown below), if G sh (l) were equal to XG sh (0) and there- 
fore G sh (n) = \ n G sh (0) for all n, C sh would be simply 
proportional to A, the adjacency matrix of the lattice 
(Aij = 1 is the vertices i and j are connected and oth- 
erwise) and the range of C sh (n) would then be limited 
to the n.n. distance. 



In order to solve the 'Ornstein-Zernike' equation (24 1, 



it is convenient to consider the simple random walk on the 
lattice where, at each time step, a particle jumps to any 
of the z neighboring sites with probability I/0. Indeed, 
G is directly related to the lattice Green function F{x) 
which is the probability generating function defined as 
(see e.g. Refs.[27||2S]) 



F ij( x ) = J2 xTp ^ 



(25) 



where pjj is the probability that the particle starting at 
i reaches j after r time steps. As is well known, one has 



F{x) = (I-*M)- 



(26) 



where M = (l/z)A. Since all sites are topologically 
equivalent, one can choose the origin of coordinates 
as the origin of the random walk and simply consider 
F n (x) = J2^Lo xT fT{ n ) where / T (n) is the probability of 
being in the nth shell after r time steps. By definition, 
Fq(x) — Fa(x) and F n (x) = c n Fij(x) where i and j are 
connected by n > 1 bonds (recall that c n = z(z— l)™ -1 ). 
It can then be shown [28] that 



F (x) = 



2(0-1) 



z - 2 + ^Jz 2 - 4(z - l)x 2 



F„ 



{x)=Cn (^l) Fq(x) iorn ^ 1 ( 2? ) 



with 



-(x) 



z - yjz 2 - 4(z - l)s 



2.r 



Using the identification 



A 



r(x) 



0-1 ' 
we readily see from Eq. (12b) that 

G sh = uF(jc) + vl 
with 

_ G sh (l) 



(28) 



(29) 



(30) 



XF (x) 
v = G sh (0) 



G 8h (l) 
A 



(31) 



Note that Eq. ( 29 ) can also be written as 

Fi(x) F (x) - 1 
xF (x) " xF Q (x) l ° ! 

which can be inverted to express i as a function of A, 

0A 



l + (z-l)A 2 



(33) 



Therefore x — > 1 when A -+ (z — 1) _1 at the spinodal. 

If v were equal to 0, one would simply have C s ' 1 = 
[uF] _1 = u _1 (I — xM) and C sh (n) would be zero for 
n > 1, as stressed above. The matrix equation C sh — 
[uF + vl]^ 1 is now easily solved: 

C sh = i[I+-F- 1 (a;)]- 1 F- 1 (a;) 
u u 

= -[i + -(i - xMjr'fi - iM) 

u u 

= — ?— [I - ^^M)]- 1 ^ - xM) 
u + v u + v 

= —^—F(x')(I-xM) 
u + v y Jy ' 



(34) 



where 



This yields 



l! 



U + V 



(35) 



C^ = -^-\F ij (x')-- r J2F ik (x')] (36) 



u + v 
where k is connected to j. Hence 

C sh {0) 



k/j 



-^ r [F (x')- X F 1 (x')}=-^-^ 

u + v z v(u + v) 



u F {x') - 1 



v(u + v) x' 



Fo(x') + - 



(37) 



and 



C sh {n) = 



1 



U + V 

1 1 



F n (x') x „ F n+1 {x') t F n _i(gQ . 

Cn Z Cn+1 Cn—1 



U + V c„ 

u F n (x') 



F n (x') - -\F n+1 {x') + (z- l)F n ^(x')} 

z 



for n > 2 



v(u + v) c„ 
where we have used the recurrence relations 



(38) 



F 2 (x') = -F 1 (x')-zF (x') 
x 

F n+1 (x r ) = -F n (x') - (0 - l)F n ^{x') for n > 2 



Introducing 



which is equivalent to 



W) - 1 

x'F (x') 



z/i 
1 + (z - l)y? ' 



(39) 



(40) 



(41) 



we finally obtain 



IV. SUMMARY AND CONCLUSION 



C sh (0) = ^F (x')(l-nx) (42) 

U + V 



and 



C sh (n) 



v(u + v) 



F (x')^ n for n > 1 . (43) 



For z = 2, one has F (x) = (1 — x 2 ) 1// ' 2 and one can 
check that Eq. ( 43 ) is in agreement with the expression 



obtained by directly solving the OZ equation in Fourier 
space. One can also check that the following susceptibil- 
ity sum-rule is satisfied, 



A § = P*"<°> + T^ 



— C"(l)]-' . (44) 



We thus see from Eq. (43 1 that C sh (n) exhibits an 
exponential decay like G sh (n), but with a different corre- 
lation length £ M = (— In |/Lt|) _1 . Moreover, since the sign 
of v = G sh (0) — G sh (l)/X depends on H, fj, is not always 
positive (see Fig. 4) and there is a range of H where the 
exponential decay is modulated by an oscillating sign. It 
turns out, however, that /i is significantly smaller than A 
in the weak-coupling (or large-disorder) regime, as illus- 
trated in Fig. 4c, so that C sh (n) decreases rapidly with 
n. Indeed, expanding all quantities in powers of J, one 
finds that 



G sh (l) = XG sh (0) + O(J 3 ) 



(45) 



so that n = 0{J 3 ). Since C sh (l) = -J/ A + 0(J 2 ), 
this implies that C sh {2) = 0(J A ) 7 C sh (3) = 0(J 7 ), etc... 
Therefore, setting C sh (n) = for n > 1 may be a rea- 
sonable approximation above the critical disorder (for in- 
stance, one has J/A c w 0.315 for z = 4 at the critical 
disorder). 

A similar calculation can be performed for the direct 
correlation function C ss (n) associated to the spin-spin 
correlation function G ss (n). It is defined via a second 
OZ equation 



*~iss f-^shf-issf-is 



(46) 



whose origin (in terms of a Legendre transform) is ex- 
plained in Ref. [10] . After some involved algebra, we 
obtain 



C S8 (n)=iJ, n - 1 [a' + b'(n-l)} 



(47) 



for n > 1, were a' , b' are functions of H / J and A/ J which 
are not detailed here for the sake of brevity. Hence C ss (n) 
has the same structure as G ss (n) with A replaced by /i 
(the fact that there is only one correlation length appear- 
ing in the final result and not two as could be expected 
from Eq. ( 46 1 is due to some remarkable cancellations 



occurring in the intermediate steps of the calculation). 
As a consequence, C ss (n) decreases rapidly with n like 
C sh (n) (i.e. C ss (2) = 0(J 4 ), C ss (3) = 0(J 7 ), etc.), so 
that the n.n. approximation is also reasonable above A c . 



In this work we have determined the two-point spin- 
spin and spin-random field correlation (or Green's) func- 
tions of the zero-temperature Gaussian RFIM on a Bcthc 
lattice along the saturation hysteresis loop. This adds 
the model to the short list of nonequilibrium systems 
for which the correlation functions are analytically cal- 
culable. In the RFIM, these functions are not known at 
equilibrium, even in one dimension, except for very spe- 
cial random-field distributions or in the universal regime 
of very small disorder. We find that the two correlation 
functions decay exponentially with the distance with the 
same correlation length. This length remains finite at 
the disorder-induced critical point, which is the expected 
behavior on a Bethe lattice. The spin-spin correlation 
function also contains a prefactor proportional to the 
distance, so that the corresponding structure factor is 
a sum of a Lorentzian and a Lorentzian squared at small 
wavevector, just like in the mean-field description of the 
equilibrium RFIM. This gives some justification for us- 
ing these simple functional forms to describe the data 
obtained in scattering experiments in RFIM-like systems 
(e.g. along the adsorption-desorption isotherms in the 
case of gases adsorbed in disordered porous solids) . We 
also find that the direct correlation functions, which are 
the inverses of the correlation functions in the sense of 
matrices, have essentially the same analytic structure as 
the correlation functions, but with a different correla- 
tion length and a modulation of the sign (depending on 
the value of the applied field). This correlation length, 
however, is small, especially in the weak-coupling/large- 
disordcr regime, and it is thus reasonable to assume that 
the range of the direct correlation functions is limited to 
the nearest-neighbor distance above the critical disorder. 
Although this 'Ornstein-Zernike' type of approximation 
breaks down in the vicinity of the critical point [55], it 
may the starting point of an analytical description of the 
hysteresis loop in the three-dimensional RFIM, as devel- 
oped recently for the soft-spin version of the model 10J . 



Appendix A: Calculation of G ss (n) and G sh (n) in the 
one-dimensional chain 



In this appendix, we present the calculation of the cor- 
relation functions G ss (n) and G sh (n) along the hystere- 
sis loop in the ID chain (specifically, along the ascending 
branch obtained by starting with a field H large and neg- 
ative). The correlations are obtained by generalizing the 
procedure used in Ref. [15] to get the magnetization. 

We first consider the spin-spin correlation function 

2 

G ss (n) = S S n — S . As stressed in the main text, 
the calculation of G ss (n) was first considered in Ref. [T3] 
but the final expression is flawed. We therefore redo the 
whole calculation, closely following the reasoning and the 
notations of Ref. [13] (note however that the calculation 
in Ref. [13] is performed along the descending branch of 



the loop). By definition, 

S ,S„ 
= $„(+, +) - $„(+, -) - $„(-, +) + $„(-, -) , 

(Al) 

where $ n (+, +) is the probability that spins at and n 
are both up, and $„(+,— ),<&„(—,+),$„(—, —) are de- 
fined analogously. 

To calculate the probabilities <&„(So, S n ) we relax the 
spins in two steps (a spin is relaxed when it is aligned 
with its local field). In the first step, the spins 5*0 and 
S n are kept down and the other spins can relax. In the 
second step, we also allow 5*0 and S n to relax. The crucial 
point is that the final state does not depend on the order 
in which spins are relaxed. 

In the first step, we need to compute the constrained 
probabilities G n (S±, S n -i) (not to be confused with the 
correlations functions) that the spins adjacent to Sq and 
S„ (see Fig. [5] for a schematic representation) are in the 
state {Si, S n -i}- 



J, Jn 






Moreover G n (— ,+) = G n (+, — ) by symmetry, and 
G n (+, +) is obtained via the sum-rule 

G„(-,-) + G„(+,-)+G n (-,+)+G, l (+,+) = 1 . (A4) 

This gives the matrix equation 

G„ = MG„_i 



where 



G n — 



G„(-,-) 
G„(+,-) 
G„(+,+) 



and 



M 



I- Pa 1 - Pi 

1 - p 1-pi 

Po 2p + pi - 1 2pi - 1 



Hence G n = M™- x Gi, with Gi 
change to the vector basis [T5] 



V 



(1-P*) 2 +1 +1 
P*(l-P*) -1 
P* 2 +1 -1 



(A5) 



(A6) 



(A7) 



We then 



(A8) 



(recall that P* = po/(l — pi + po)) in which the matrix 
M takes the form 



FIG. 5: Schematic representation of the environment of the 
spins So and S„- In addition to the spin variables Si = ±1, 
we also use the variables h — (1 + S-i)/2, li — (1 + <Si)/2, 
r2 = (1 + S n _i)/2, and ri = (1 + S' n +i)/2 taking the values 
0,1. 



M = V^MV = 



1 

pi - po 1 - Pi 
Pi-pa 



so that 



(A9) 



1. Constrained probabilities G„(Si,S„-i) 



By definition, 

G„(<Si, Jn-l) = 



S 2 ...S„_: 



P(5i,S 2 ...5 n _ 2 ,5 n _i) (A2) 



where P(Si, S2 ■ ■ ■ 5'n-2i SVi-i) is the probability of the 
configuration {Si, S*2 . . . 5 f n _2, Sn^i} when the spins So 
and S„ are pinned down and the spins between them are 
allowed to relax. 

Carrying out the calculation as in Ref . [T3] , one eas- 
ily derives the following recurrence relations for the con- 
strained probabilities along the ascending branch of the 
hysteresis loop, 

G„(- -) = (1 - p )G n _i(- -) + (1 - pi)G n _i(- +) 

G„(+, -) = (1 - p )G n _!(+, -) + (1 - pi)G„_i(+ J +) . 

(A3) 



M 



1 

" _1 = (pi-poT- 1 (n - 1)(1 - Pl )(pi - Po ) n - 2 

(pi-poY 1 - 1 

( A1 °) 

To simplify the notation, we shift to the variables I2 = 
(1 + Si)/2 and r 2 = (1 + S„_i)/2 that take the values 
0, 1, and after some algebra we finally obtain 



G„(0,0) = (1-P*) 2 +[P*(1-P*) 

1 - -Pi , n* 



(n- 1)P* 



P*](pi-Po) r 



Pi -Po 
G„(1,0) = P*(1-P*)-[P*(1-P*) 

+ («-i)p*^^ L ]bi-p r- 1 

Pi -Po 

G„(1,1) = P* 2 +[P*(1-P*) 

+ („ _ l)p* I^IL _ P *l ( pi _ po) n-l 



Pi -Po 



(All) 
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2. Calculation of G 3S {n) 

To compute <&„(So, S n ) and then G 3S (n) we now con- 
sider the second step where So and S n are relaxed. We 
define P(So, S n \l±, Si . . . S n -i, r±) as the probability of 
the state {So,S ra } under the condition that the state 



{li,Si . . . S n -i,ri} has been reached after the first step 
(as indicated in Fig.[5| h and n describe the states of the 
spins S_i and S n +i, respectively). Knowing these prob- 
abilities and the probability of all possible environments 
of So and S n after the first relaxation step, we can write 



$ n (S ,S n ) = J2 P * h+ri ( 1 - p *f~ h ~ ri E P(Si---S n -i) Yl P(S 1 ...S n -. 1 )P(S ,S n \l 1 ,S 1 ...S n - 1 ,r 1 ). 

Jin Si...S„_i Si...S„_i 

(A12) 



r 



appear which account for the cases where the flip of So 
(resp. S n ) triggers an avalanche which make all the spins 



If any of the spins {Si, . . . , S n -i} is up the probabil- 
ity related to the second relaxation step is just a prod- 
uct of two independent terms. On the other hand, if {Si...S n _i} to flip up, changing the environment or the 
all the spins {Si, ■ ■ . , S„_i} are down, the expressions of state of S n (resp. Sq). This yields 
the probabilities are more complicated since extra terms 



P(-l,-l|ii,5i ...£„_i,ri) = (1 - p h+ i 2 )(l - p ri+r2 ) 



P(+l J -l|li > Si...S B _i,ri) 



P(-l,+l|Ji,5i...S B _i,n) = 



P(+l,+l|J 1 ,Si...S„_i,n) 



Ph+h( l -Pr 1 +r 2 ) 

( - \ n ~ 1 
Phi 1 ~ PrJ ~ [T=%r ) PhiPn+l ~ Prx) 

(1 - p h+ i 2 )p ri+r2 

(l-p h )p ri - {tj^J (Ph+1 ~Ph)Pry 

Ph+l 2 Pri+r 2 



Pi 



iPn + (tt^ 1 ) [Ph(p ri +i -Pn) + (Ph+i -Ph)Pn] H {Si . . . S n -i} = {-1 



if {Si . 


•S„_i}?£{-1-- 


-1} 


if {Si . 


.S n _i} = {-1-- 


-1} 


if {Si . . 


.S n _i}^{-1--- 


-1} 


if {Si . . 


.S„_i} = {-1--- 


-1} 


if {Si . 


..S n _i}^{-1-- 


•-1} 


if {Si . 


..S n -i} = {-!■■ 


•-1} 



where l 2 = (1 + S x )/2 and r 2 = (1 + S n _ x )/2. Using 
Eqs. pS| , ( |A13| ), and the probability (l-p )"~ 1 that all 
the spins {Si . . . S n _i} are down after the first relaxation 



(A13) 



I 

step, we then obtain 

*«(-.-) = E p (^) 

h,ri 



(2,r 2 

*n(+,-) = YP{li,ri){^2G n {h,r2)pi 1+ i 2 {l-p ri +r 2 ) 
h,ri l 2 ,r 2 

- {Pi ~ PoT~ l Ph(Pr x +l -Pn)\ 

*n(~, +) = X] P (^' r l){ 53 G »^ 2 ' r2 )( X _ Pli+b )Pr! 
Ji.ri l 2 ,r 2 

- (Pi -Po) n-1 (pji+l -P/i)Pn] 

*«(+»+) = 53 p (^: r i){5Z Gn ( /2 ' r2 )^i+ i ^'-i 



+ '('2 



I_+I"2 



\n-l„ 



+ (Pl-Po)™ PZi (Pn+l -Pn) 

+ (Pi -Po)" _1 (P/i+l -PzJPn}- 



(A14) 
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where P (l, r) = P* l+r 

the following exact relations are satisfied 



(1-P*) 2 ' r . One can check that 



*„(+,+) +2* n (+,-) 

*n(+,+) 



*n(+,-)=*n(-+) 
+ $„(-,-) =1 

+ *„(+, -)=i(m(ff) + l) 



(A15) 



where the magnetization m(H) is given by Eq. Q. 

In Fig. pi Eqs. (A14) are compared to simulation re- 
sults in the case n = 4. The excellent agreement confirms 
that the whole calculation is correct. 
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FIG. 6: Probabilities $ n (So, S4) along the ascending branch 
of the hysteresis loop for z = 2 and A = 4. The simulation re- 



sults (symbols) are compared to the predictions of Eqs. ( A14 1 
(lines). Simulations were performed on random graphs with 
N = 10 6 and the results were averaged over 1000 disorder 
realizations. 



A simpler expression is actually obtained by using Eqs. 
( |A15 | to express G ss (n) in terms of $„(—,—) only, as 
done in Ref.|13|. This finally yields 



G l 



*(n) = 4{*„(-,-) 



.1 -m(H) 



( Pl ~Po) n -"[a + b{n-l)] (A17) 



with 



a = 4P*{Q* - P*) [2(1 - P*) - P*(Q* - P*)} 
^2 (1-Pi) 



b = 4P*{Q* -p*y 



For z — 2 we recall that 



P* = 



(pi - Po) 



(A18) 





Po 


1- 


- Pi + Po 


Pl 


- Pi + P0P2 



1 - Pi + Po 



(A19) 



3. Calculation of G sh (n) = S K 



We now consider the spin-random field correlation 
function 



Note that in Ref.[T3] it is stated that the proba- 
bilities $> n (So,S n ) are just linear combinations of the 
constrained probabilities G„'s, without inhomogeneous 
terms. Eqs. 

*n( — ) 



(A14) show that this is only true for 



Finally, inserting Eqs. (|A14|) in Eq. ( |Al| ) yields 

2pz 1+ * 2 )(l 



S S n = J2 P(h,r x ) J2 Gn(h,r 2 ){l 

h,ri h,r 2 

-1 V^ 



G 8h {n) = S h r , 



E 

So,S„ 



dh n Soh n Q n (So, S n , h n ) 



2p ri +r 2 ) 



+ 4(pi - poY 



h,ri 



P{h,ri)p h {p r 



■Pri) ■ 



(A16) 



(A20) 

where $„(So, S n , h n ) is the probability density that the 
two spins So and S n are in the state {So,S n } after full 
relaxation of the system, with the random field acting on 
the spin S n having a value within (h n ,h n + dh n ). The 
calculation of these probabilities is straightforward since 
we already have computed the probabilities $„(So, S n ). 
We only need to restrict the integration of the random 
field distribution p(h) over a range of h compatible with 
the state of the spin S n . This gives 
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(-, -, h n ) = J2 P(h,ri){ Y, G «^ 2 ' ^H 1 - Ph+i 2 )p(h n )Q[h n < 2(1 - ri - r 2 )J - H}} 

h,ri h,r 2 

(+, -, h n ) = J2 P(h,n){ J2 G n(h, r 2 )p h+ i 2 p(h n )@[h„ < 2(1 - n - r 2 )J - H] 

- (pi - p ) n_1 P/ 1 /3(^)e[2(l _ n - l)J - if < ft n < 2(1 - n)J - H]j 
*«(-, +, M = ^ P(ii,ri){ J] G n(^, r 2 )(l - p, 1+Ia )p(/i„)6[/i n > 2(1 - n - r a ) J - if] 

- (Pi -Po) n_1 (p/ 1+ i - PJx)p(An)»[/in > 2(1 - r x )J - H}} 

*»(+, +, /in) = J! P (^l){ X! G "(' 2 ' r 2)pi 1+ i 2 p(ft„)e[/l n > 2(1 - n - r 2 )J- H] 



h,r\ l 2 ,r 2 

\n-l„ 



+ (pi -p )""V/ 1 p(^n)e[2(i - n - i)J - F < K < 2(1 - n)j - H) 

+ {pi-p ) n -\pi 1+ i-Ph)p(h n )@[h n > 2(1 -n)J- H}} 

I 



where 9(.) is the characteristic function of the domain outside). Inserting in Eq.(A20l yields 
indicated by the argument (i.e. 1 inside the domain and 



(A21) 



S a h n = / dh n K[$n(+,-,K) + ®n(+,+,h n ) - $n(-,-,h n ) - $ n (-,+,hnj\ 

/+oo 
dh n h n p(h n ) 

«i,n l2 ,r 3 -°° 

r+oo 

+ 2(pi -po)™ -1 V" P(ii,ri)(pj 1+ i -PiJ / dh n h n p(h n ), 

I. ... J2(l- ri )J-H 



(A22) 



and finally 



S'o/in = 2A(pi -p ) r 



^P*^(l-P*)i-'i(p, 1+1 



Pi,, 



where Q* is given by 

°" = E f z I x ) i^wf [i - p*(if)]^ 1_fc p fc+ i(if) 



fe=0 



^ p*n (1 _ p*)l-r lp ( 2(1 _ ri) J _ jgr) 



and 



(B2) 



(A23) 



fe=0 



G sh (0) = 2A£ (*) [P*(H)] k [1-P*(H)Y 



Appendix B: Calculation of G ss (2) and G Bh {l) on the 
Bethe lattice. 



p((z-2k)J-H) 



(B3) 



We recall that P* (resp. Q*) is the probability that, 
along the ascending branch of the loop, a spin is up given 
In this Appendix we calculate G ss {2) and G sh {l) on a that a neighbor is forced to be down (resp. up). 



Bethe lattice with coordination number z. For complete- 
ness, we first recall the expressions of G ss (l) and G sh (0) 
obtained in Ref.jHj: 



1. Calculation of G ss (2) 



G ss {l)+m 2 = 1-4P* +AP* 



To compute the correlations between two spins S\ and 
(Bl) S2 at the distance n — 2 we consider a central spin So 
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and its z neighbors {Si, . . . S z } as depicted in Fig. [7^,. 
By definition, 

S1S2 — 2_^ SiS 2 p(So, Si . . . s z ) 

Sq,Si...S z 

= E S ^i E P(-l,Si...S z ) 

Si,S2 Sz---S z 

+ J2 P(+l,Si...S z j\ 

S 3 -S z 
= Y, SiS 2 [P(-l,Si,S 2 ) + P(+l,Si,S 2 )} (B4) 

Si,S2 

where P(So, 6\, . . . , S z ) is the probability of having the 
configuration {So, Si, . . . ,S Z } when the system is fully 
relaxed. 



forced to be down). In Eq. (B6), the summation ac 



counts for all the different ways of having Sq = 1 when 
q neighbors are up, and the term in front of the summa- 
tion accounts for the probability that So has z — q of his 
neighbors down if it is up. 



Using Eqs. (B5) and (B6) we find 



P(— 1, Si, S 2 



^( z , 2 )[i-pr'-"[Ff"(i- ;)i;i! 



1=0 

z-2 



1=0 



I 

z-2 



l+n 



p( + i,5 llS2 )=E(^i[i-Qr-[Ef + ; 



fe=0 



where 



[Q* - P*] l+n ~ k [P*] k Pk ] 



Si+S 2 + 2 
n = = 0, 1, 2 



(B8) 



(B9) 




We finally obtain the following expression for the corre- 
lations at the next-nearest neighbor distance 



^=E(-W!)e V: 2 



71=0 



nl ' — ' \ I 
1=0 



x [i-p*r'- n in l+n (i-Pi + n) 



(/:,) 




z-1 



z-1 



FIG. 7: Schematic representation of the environment of (a) a 
single spin So and (b) a pair of spins Si and Sj . 



By relaxing the spins in two steps, we obtain 

p(-i, Si...s z ) = [i- p*r q m q (1 - *>,) (B5) 



p(+i,si...s z ) = [i-Q*r q J2 



k=0 



[Q* - P*] q [P*] Pk 



(B6) 



where q is the number of neighbors of So that are up 
when the system is fully relaxed, 



ELi Sj + z 



= 0,l,...z. 



(B7) 



Eq. ( |B5[ ) is rather straightforward (recall that P* is the 
probability that a spin is up given that a neighbor is 



[1-01 



1-1-7 



l+n 



Y,( l t n )iQ*-p*} l+n - k [p*] k p* 



Lfc=o 



k 



(BIO) 



2. Calculation of G 3h (l) 



To compute 
G sh (l) = S~h] = J2 f dhjSihjPiSi, hj) , (Bll) 

we again relax the spins in two steps so to obtain 

p(s t ,h,) = j2( z ] 1 )p* l (i- n*- 1 - 1 
1=0 ^ ' 

^J2( z ' r 1 )p* r a-n z - 1 - r p(s l ,h J \i,r), 

(B12) 

where P(Si, hj\l, r) is the probability that the spin at i is 
in the state Si after the second relaxation step, with the 
random field acting on Sj having a value within (hj, hj + 
dhj), and under the condition that the environment of 
Si and Sj is in the state (I, r) after the first relaxation 
step (see Fig. [7]b). As in Ref.[8j, I = 1, . . . , z — 1 (resp. 
r = 1, . . . , z — 1) is the number of neighbors of Si (Sj) 
that are up, without taking into account Sj (resp. Si). 
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These probabilities are given by 
P(-l,h 3 \l,r) : 
P(+l,h 3 \l,r)-- 



For z — 2 and n = 1, one can check that this equation 



(1 - pi)p{hj) if hj <(z- 2r)J - H 
(1 - pi+i)p(hj) if /ij > (z - 2r) J - # 

Pip(hj) if /ij < (z — 2r)J — H 

pi + ip(hj) if frj > (z - 2r) J - if . 

(B13) 



gives back Eq. (A23) 



Inserting these expressions in Eq. ( B12 1 and after some 
algebra we finally obtain 



G sn (l) = 2A 



E(V)^(i- /,i) ' 'U',^-i>n 



2=0 



I 



E 



z- 1 



p* r /1 p* 



r p{{z-2r)J -H) 
(B14) 
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